## MODELING ENDOGENOUS EXPLANATORY VARIABLES

model {
	for (i in 1:N) {                     # N OBSERVATIONS
	    y[i] ~ dnorm (yhat[i], tau.y)
	    yhat[i] <- inprod(B[], X.1[i,]) + D[country[i],1]   
	    
	    w1[i] ~ dnorm (w1.hat[i], tau.w1)
	    w1.hat[i] <- inprod(G.1[], W.1[i,]) + D[country[i],2]
	    
	    w2[i] ~ dnorm (w2.hat[i], tau.w2)
	    w2.hat[i] <- inprod(G.2[], W.2[i,]) + D[country[i],3]	    

	    w3[i] ~ dnorm (w3.hat[i], tau.w3)
	    w3.hat[i] <- inprod(G.3[], W.3[i,]) + D[country[i],4]	  

	    w4[i] ~ dnorm (w4.hat[i], tau.w4)
	    w4.hat[i] <- inprod(G.4[], W.4[i,]) + D[country[i],5]
	    
	    w5[i] ~ dnorm (w5.hat[i], tau.w5)
	    w5.hat[i] <- inprod(G.5[], W.5[i,]) + D[country[i],6]	    

	    w6[i] ~ dnorm (w6.hat[i], tau.w6)
	    w6.hat[i] <- inprod(G.6[], W.6[i,]) + D[country[i],7]	
    }

	for (j in 1:J) {
		D[j,1:K] ~ dmnorm (mu.D[1:K], Tau.D[1:K,1:K])
	}

	for (k in 1:K) {
		mu.D[k] <- 0 
	}
	
	# COVARIANCE MATRIX IN JOINT DISTRIBUTION
	Tau.D[1:K,1:K] ~ dwish (W[1:K,1:K], K+1)
	
	Sigma.D[1:K,1:K] <- inverse(Tau.D[1:K,1:K])

	# COMPUTE THE VARIANCE TERMS	
	sigma2.d <- Sigma.D[1,1]
	sigma.d <- sqrt(sigma2.d)

	sigma2.22 <- Sigma.D[2,2]+sigma.w1^2
	sigma.22 <- sqrt(Sigma.D[2,2])

	sigma2.33 <- Sigma.D[3,3]+sigma.w2^2
	sigma.33 <- sqrt(Sigma.D[3,3])

	sigma2.44 <- Sigma.D[4,4]+sigma.w3^2
	sigma.44 <- sqrt(Sigma.D[4,4])

	sigma2.55 <- Sigma.D[5,5]+sigma.w4^2
	sigma.55 <- sqrt(Sigma.D[5,5])

	sigma2.66 <- Sigma.D[6,6]+sigma.w5^2
	sigma.66 <- sqrt(Sigma.D[6,6])

	sigma2.77 <- Sigma.D[7,7]+sigma.w6^2
	sigma.77 <- sqrt(Sigma.D[7,7])

	# DERIVE THE COVARIANCE			
	cov.12 <- Sigma.D[1,2]	
	cov.13 <- Sigma.D[1,3]
	cov.14 <- Sigma.D[1,4]	
	cov.15 <- Sigma.D[1,5]	
	cov.16 <- Sigma.D[1,6]
	cov.17 <- Sigma.D[1,7]	
	
	# COMPUTE CORRELATION
	rho.12 <- cov.12/sqrt(sigma2.d*sigma2.22)
	rho.13 <- cov.13/sqrt(sigma2.d*sigma2.33)
	rho.14 <- cov.14/sqrt(sigma2.d*sigma2.44)
	rho.15 <- cov.15/sqrt(sigma2.d*sigma2.55)
	rho.16 <- cov.16/sqrt(sigma2.d*sigma2.66)
	rho.17 <- cov.17/sqrt(sigma2.d*sigma2.77)

    # PRIOR OF THE PRECISION PARAMETER
    tau.y ~ dgamma (0.01, 0.01)   # dgamma (r, mu), WHERE r IS THE SHAPE AND mu IS THE RATE PARAMETER      
    sigma.y <- 1/sqrt(tau.y)  

   	tau.w1 ~ dgamma (0.01, 0.01)      
    sigma.w1 <- 1/sqrt(tau.w1)   

   	tau.w2 ~ dgamma (0.01, 0.01)      
    sigma.w2 <- 1/sqrt(tau.w2) 

    tau.w3 ~ dgamma (0.01, 0.01)      
    sigma.w3 <- 1/sqrt(tau.w3) 

   	tau.w4 ~ dgamma (0.01, 0.01)      
    sigma.w4 <- 1/sqrt(tau.w4)   

   	tau.w5 ~ dgamma (0.01, 0.01)      
    sigma.w5 <- 1/sqrt(tau.w5) 

    tau.w6 ~ dgamma (0.01, 0.01)      
    sigma.w6 <- 1/sqrt(tau.w6) 
    
    # PRIORS ON UNMODELED COEFFICIENTS
    for (k in 1:M.1) {
    	B[k] ~ dnorm (0, 0.0001)
    }

    for (l in 1:L.1) {
    	G.1[l] ~ dnorm (0, 0.0001)  
    }   

    for (l in 1:L.2) {
    	G.2[l] ~ dnorm (0, 0.0001)  
    }   

    for (l in 1:L.3) {
    	G.3[l] ~ dnorm (0, 0.0001)  
    }   

    for (l in 1:L.4) {
    	G.4[l] ~ dnorm (0, 0.0001)  
    }   

    for (l in 1:L.5) {
    	G.5[l] ~ dnorm (0, 0.0001)  
    }   

    for (l in 1:L.6) {
    	G.6[l] ~ dnorm (0, 0.0001)  
    } 
       
} # END OF MODEL

